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Parafermions are exotic quasiparticles with non-Abelian fractional statistics that can be realized and stabilized 
in 1-dimensional models that are generalizations of the Kitaev p-wave wire. We study the simplest generaliza¬ 
tion, i.e. the Z 3 parafermionic chain. Using a Jordan-Wigner transform we focus on the equivalent three-state 
chiral clock model, and study its rich phase diagram using the density matrix renormalization group technique. 

We perform our analyses using quantum entanglement diagnostics which allow us to determine phase bound¬ 
aries, and the nature of the phase transitions. In particular, we study the transition between the topological and 
trivial phases, as well as to an intervening incommensurate phase which appears in a wide region of the phase 
diagram. The phase diagram is predicted to contain a Lifshitz type transition which we confirm using entan¬ 
glement measures. We also attempt to locate and characterize a putative tricritical point in the phase diagram 
where the three above mentioned phases meet at a single point. 


Introduction — There has been concerted effort to engineer 
systems with stable Majorana bound states, and other anyonic 
quasiparticles, for use in the topological quantum computa¬ 
tion architecture [1-7]. For example, there has been recent 
progress in attempts to isolate Majorana bound states in quan¬ 
tum nanowires [5, 8-10] and in superconductor surfaces im¬ 
planted with a line of magnetic impurities [11]. These quasi- 
ID systems effectively realize a version of the Kitaev p-wave 
wire model [12], and are predicted to have a gapped topologi¬ 
cal phase which supports characteristic Majorana bound states 
at the ends of the wire. 

While the boundary modes in these heterostructure systems 
are non-Abelian anyons, they are unfortunately known to be 
insufficient for universal quantum computation. A possible 
remedy for this problem has been to look for more exotic non- 
Abelian excitations. For example, Fendley has recently sug¬ 
gested exploring one-dimensional para-fermionic models 
which support topological phases with more computationally 
efficient non-Abelian any on bound states [13]. Still, the Zn 
non-Abelian anyons are not able to perform universal quan¬ 
tum computation, however they can be leveraged to create a 
2D phase with Fibonaccci anyons, which are universal [14]. 
These promising features have spurred wide spread interest 
in these models, and has led to many analytical and numerical 
studies, including several experimental proposals for realizing 
these topological phases [15-39]. 

In this work, we continue along these lines of research by 
exploring the rich phase diagram of the Z3 para-fermionic 
chain; though for ease of calculation we actually study the 
Jordan-Wigner transformed para-fermionic chain[40], includ¬ 
ing chiral interactions. The resulting model is the three state 
chiral clock model. This model re-surfaced in this context in 
Ref. 13 as a candidate for exhibiting non-Abelian bound states 
beyond Majorana fermions. It was shown analytically that 
para-fermionic boundary zero modes can exist in this model 
when spatial-parity and time-reversal symmetries are broken 
via chiral interactions [13]. This was verified numerically in 
Ref. 41, which confirms that chiral interactions can help to 
stabilize the boundary zero modes, although the zero modes 
themselves are more fragile than one might initially expect. 


Here we are interested in studying the full phase diagram of 
the chiral clock model as a function of two chiral-interaction 
phase-parameters (#, </>), as well as the relative strength of the 
nearest neighbor coupling ( J) to the local Zeeman field (/). 
Using entanglement techniques, we have been able to locate 
the phase boundaries that separate the topological phase from 
the trivial gapped phase, and a critical incommensurate phase, 
the latter of which has no analog in the Kitaev p-wave wire 
model. We have conclusively identified the region in which 
there is a topological phase, and have explored the nature of 
the quantum phase transitions in and out of the three adjoining 
phases. In addition, by studying oscillatory properties of the 
system in, or near, the incommensurate phase, we establish the 
approximate location of a putative tricritical point[42, 43], and 
further support the entanglement signatures that were recently 
proposed for identifying Lifshitz transitions [44]. 

The article is arranged as follows. We first discuss the de¬ 
tails of the model, and the criteria used to map out its phase 
diagram. For our numerical simulations, the density matrix 
renormalization group (DMRG) [45, 46] algorithm is em¬ 
ployed, as it gives immediate access to the entanglement en¬ 
tropy (EE), and therefore the central charge, at putative critical 
points/regions in the phase diagram [47] . Next, we discuss the 
general features of the phase diagram and locate regions in the 
topological phase (where para-fermion boundary modes may 
exist). We also discuss the nature of the phase transitions out 
of the topological phase. For part of our study we discuss our 
observations pertaining to a critical incommensurate phase, 
and the possibility of a tricritical point [42, 43] in the phase 
diagram at the intersection of the topological, trivial, and in¬ 
commensurate phases. We also find a region of the phase di¬ 
agram which exhibits the critical entanglement features of a 
Lifshitz transition [44]. Finally, we conclude by discussing 
future directions and possible relevance to experiments look¬ 
ing for para-fermions. We also include four appendices which 
discuss some subtleties of the numerical analysis. 

The Model — For our study we use the ID 3-state (Zf) chiral 
clock model [13, 42, 48, 49]. The Hamiltonian for the 3-state 
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chiral clock model is: 

L L -1 

Hs = -fJ2 r}e~^ a } a j+ ie ~ ie + h -c- (1) 

3 = 1 3 = 1 

following the notation in previous work [13], where /, J, 6 
and (j) are scalar parameters, and ai and r- L are local three state 
spin operators on site i. The spin operators have the properties 
7~ 3 = a 3 = /, err = u) t(j , where cc = e 2 W 3 # Specifically, 
we use the matrix representation 

/ 1 0 0 \ /0 1 0 \ 

T = 0 U3 0 , < 7=1 0 0 1 J . 

\0 0 U ) 2 J \1 00 / 

The chiral clock model is related to the para-fermionic 
chain through a Jordan-Wigner transformation [13, 40], simi¬ 
lar to the well-known, analogous case that the Kitaev p-wave 
wire is related to the transverse-field Ising model via the same 
type of transformation. The parafermion operators are defined 
as 



at site j. The corresponding para-fermionic Hamiltonian is 

L L -1 

H 3 = -fJ2 - Ju 1 Y, V’jxi+ie-" + h.c. (4) 

3 =1 3=1 

The chiral clock model has a global Z 3 symmetry that can 
be represented with x = Ylf=i 7 j = Here Z is the 

generator of the symmetry, and has three different eigenvalues 
0,1,2. In addition, when all of the coefficients of the Hamilto¬ 
nian are real, i.e. when the system is a Z 3 -ferromagnet or Z 3 - 
anti-ferromagnet Hamiltonian, then the Hamiltonian is invari¬ 
ant under time-reversal, charge-conjugation, and parity sym¬ 
metries. This can be easily seen from the following definitions 
of these symmetries. Charge conjugation C acts on the spin 
operators via C<jjC = cc 2 crj, CrjC = rj, C 2 = 1 . As an 
aside, note that charge conjugation, together with the Z 3 sym¬ 
metry, forms the S 3 permutation symmetry, i.e. the symmetry 
obeyed when the 3-state clock model is restricted to the 3- 
state Potts model. Time reversal T acts on the spin operators 
via TdjT = Gj , TtjT = rj, T 2 = 1, and complex conju¬ 
gates any scalar coefficients. Spatial parity P acts on the spin 
operators via PgjP = cr_j, PrjP = r_j, P 2 = 1. Finally, 
we note two things: (i) due to the symmetry of the Hamil¬ 
tonian with respect to (j) and 6 , we only need to consider the 
region of the phase diagram where </> and 0 each range from 0 
to |, and (ii) for / = J, the system is self-dual along the line 
<p = 0. The details of these two properties are in Appendix A. 


There are many previously known results about this model 
(Eq. 1), beginning with the original proposals of Ostlund [42] 
and Huse [48]. For example, the corresponding two- 
dimensional classical Hamiltonian for 0 = 0 was studied in 
Ref. 42, and the the one-dimensional quantum Hamiltonian 
was studied in Ref. 43 for the restricted case i> = 0. One of 
the most important early results is that Eq. (1) has a second 
order quantum phase transition at / = J when 0 = 4> = 0. At 
this point the model realizes the full S 3 permutation symme¬ 
try (instead of just Z 3 ), and the critical point is described by 
the critical conformal field theory for the 3-state Potts model, 
which has central charge 4/5 [50]. In addition, the line 
/cos(3</>) = Jcos(3$) [51] is known to be integrable and 
4> = 0 = | is super integrable [52, 53]. Despite this, the 
knowledge of the location of some important critical points 
and their associated properties is an open question. 

Generically, it is known that the phase diagram is divided 
up into two gapped regions, one of which is identified with 
small values of / (compared with J), and the other with large 
values of /. These regions are separated by continuous quan¬ 
tum phase transitions that we will identify and discuss further 
below. Using a more modem terminology, the gapped phase 
for small / is a symmetry broken phase of the 3-state clock 
model and it exactly corresponds to the “topological" phase 
in the Jordan-Wigner transformed para-fermionic chain. The 
gapped phase for large / is a disordered phase of the 3-state 
clock model, and maps onto the “trivial" phase of the para- 
fermionic chain. This gives another example of a case where 
the degeneracy associated to symmetry breaking is mapped 
to topological degeneracy via the Jordan-Wigner transforma¬ 
tion [54, 55]. Hence, in either representation this phase has a 
three-fold ground-state degeneracy, which can be detected by 
measuring the ground state EE. On the other hand, the triv¬ 
ial phase is equivalent to the spin disordered phase, which 
does not have a generic ground-state degeneracy. The param¬ 
eter / is thus an important tuning parameter for the phase di¬ 
agram, and analogous to the external transverse field in the 
Ising model. 

While we expect these general features to pervade the phase 
diagram, the phase space for generic 0 and 0 is largely un¬ 
explored. Additionally, it is known that the combination of 
the Z 3 , symmetry and the chiral nature of the interactions, 
gives rise to interesting behavior that cannot be found in the 
Majorana/Ising case. For example, this model supports a so- 
called “incommensurate phase" which is not present in the 
transverse-field Ising model with chiral interactions [42]. 

This motivates the main objective of our article, which is 
to characterize the phases and the nature of the phase tran¬ 
sitions over the entire phase space. We will show that there 
are two types of phase transitions that occur to destabilize the 
topological phase, and there is a large region of critical incom¬ 
mensurate phase that separates the topological from the trivial 
phase over a wide range of parameters. Let us now move on 
to a discussion of the methods we employ. 

Methods — We primarily use the spatial EE in order to char¬ 
acterize the phase diagram. This measure has been widely 
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used to detect topological order in 2D [56, 57], and has been 
applied more recently to ID topological phases [58]. The EE 
can be derived by partitioning the system into two regions A 
and T?, and then calculating the reduced density matrix of re¬ 
gion A by tracing over all the degrees of freedom in region 
B. Mathematically, the reduced density matrix is given by 
pA = Tr Bp, and the corresponding entanglement entropy is 
defined to be: 

S = -Tr (p A In pa) • (5) 

There are two useful entanglement indicators we will em¬ 
ploy to identify the phases and phase transitions for the chiral 
clock model. First, for the gapped regions of the phase di¬ 
gram, it is known that for one dimensional gapped systems 
the entanglement entropy increases with the the block size 
l (the size of region A), and saturates when l reaches the 
correlation length [47]. Furthermore, if there is topological 
ground-state degeneracy we would expect an entanglement of 
order ~ logD where D is the degeneracy [58]. To eliminate 
the most harmful finite-size effects we will take the central- 
cut, i.e., cutting the chain in half, to identify the nature of the 
gapped phases. 

For critical regions of the phase diagram, it is known 
that the entanglement entropy will grow logarithmically with 
system size, and the scaling is characterized by the central 
charge [47]. More specifically, for critical systems with open 
boundary conditions, the form of the entanglement scaling law 
is [47]: 

c f2L . 7rl\ 

s = 6 ln v V sm T) + s ° 6 

where l is the length of the subsystem, c is the central charge, 
and So contains the sub-leading corrections. Once we know 
the central charge we will have an important piece of infor¬ 
mation about the phase transition/critical phase, and can then 
appeal to previously known analytic results in restricted parts 
of the phase diagram to help further specify the phase dia¬ 
gram. Below we will see the efficacy of these two indicators 
for determining the phase diagram. 

To arrive at our results for the phase diagram (and to ob¬ 
tain reasonable estimates of the phase boundaries in the ther¬ 
modynamic limit), we simulated Hamiltonians using open¬ 
boundary DMRG with 100 sites, and a bond dimension m = 
100. We find this to be sufficient for the phases with low en¬ 
tanglement entropy. For the critical phases, additional checks 
were performed with bond dimension m = 200. For estab¬ 
lishing characteristics of other phases, for example, the region 
of critical incommensurate phase, larger lengths of 400 sites 
were also tested. 

Results — Fet us now move on to discuss the results of 
our numerical calculations. First, we present the full three- 
parameter phase diagram (/,0,0) over the reduced domain in 
Fig. 1, where we have set J = 1—/. The basic topology of the 
phase structure is clear. We find three distinct phases as men¬ 
tioned above. The phase corresponding to largest / values 



Figure 1. Three-dimensional phase diagram of the chiral 3-state 
clock model in terms of /, 6 and 0 with J — 1 — /. For details of 
the Hamiltonian see Eq. (1). The topological, trivial, and incommen¬ 
surate (IC) phases are indicated. The coloring is a function of the 
value of / at the critical surface separating the phases. The dashed 
line that connects points (0, 0, 0.5) and (7t/3, 7t/3, 0.5) is the self¬ 
dual line. 


is generically the trivial phase, and the phase corresponding 
to the smallest / values is generically the topological phase. 
They share a common/direct phase boundary between them 
when 6 and 0 are small. For large 6 or 0 , an intermediate 
incommensurate phase appears between the two. 

We show the central-cut EE in Fig. 2(a),(b),(c) for several 
2D cross-sections of the 3D phase diagram. These plots help 
to identify the gapped phases and the topology of the phase 
boundaries. To more clearly identify the nature of the critical 
regions/boundaries we also calculate the central charge via the 
scaling relation. It is interesting to see that the observed loca¬ 
tions of the phase boundaries for cross sections 0 = 0 and 
0 = 0 are broadly consistent with earlier works [42, 43], and 
that the topological phase itself is stable over a large part of 
the phase diagram [60]. 

We indicate several special points on these cross sections: 
Point A in Fig. 2(a) and Fig. 2(c) is the transition point of 
the three-state Potts model associated with c = 4/5 [50], 
and Point B and C are putative tri-critical points. We indi¬ 
cate approximate locations of the phase boundaries with solid, 
dashed, or dot-dashed lines, depending on the nature of the 
phase transition, as indicated in the figure caption. Finite size 
effects were checked around specific points along the criti¬ 
cal lines by running system sizes of L -100 to 400 on a finer 
grid. The locations of these lines did not change significantly 
in comparison to the resolution of our grid, except in certain 
regions which are discussed in further detail in later sections. 

From the central-cut EE we see that the trivial phase is char¬ 
acterized by a small EE, while the topological phase has a 
nearly uniform EE of ~ ln3 indicating a three-fold degen¬ 
eracy of the ground state. The change of EE is abrupt be¬ 
tween the two phases as can clearly been seen in Fig. 2(a) and 
Fig. 2(c) for 6 < 7r/4 and 0 < tt/ 6 respectively. We also 
verified that this transition is accompanied by a divergence 
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Figure 2. Three cross-sections corresponding to (a) 4> = 0 (b) 0 = 7t/3 and (c) 4> — 0 of the three dimensional phase diagram, and all for 
L — 100. Topological, trivial, and incommensurate (IC) phases are identified by the central-cut entanglement entropy (color coded). For (a) 
and (b) a 2D grid in increments of 0.01 was used to resolve fine features of the transitions, (c) was mapped out on a 2D grid in increments of 
0.05. Point A is the transition point of the 3-state Potts model, i.e. the chiral clock model for (0 = (j) — 0). Points B and C are Lifshitz points 
and are associated with putative tricritical behavior. The solid lines, dashed lines, and dotted-dashed lines indicate direct topological-trivial 
(c = 4/5) type, Kosterlitz-Thouless type, and Pokrovskii-Talapov type [59] transitions respectively. The thick circularly-dotted line represents 
an upper bound on the region where exact parafermionic zero modes can exist [13]. Panels (d), (e) and (f) show the corresponding central 
charges for cross sections (a),(b),(c) respectively. The IC phase is associated with central charge c — 1 (yellow) whereas the critical regions 
close to point A have c = 4/5 (green). 


in the second order derivative of the ground state energy (not 
shown). 

The third phase in the phase diagram is the incommensurate 
phase. This is a critical phase in which the correlation func¬ 
tions generically behave as A(r)e^ 27r2 / 3 ^ r , where A decays 
algebraically and Q is irrational. The oscillatory properties 
of the correlation functions also manifest themselves in oscil¬ 
latory behavior seen in energy gaps, which we address later. 
Although there is not an extremely sharp distinction between 
the central-cut EE for the topological and incommensurate 
phases, the EE scaling with system size is markedly differ¬ 
ent. The former has an EE that quickly saturates to a constant 
value of In 3 with sub-system size, while the latter has EE that 
diverges logarithmically with sub-system size. By fitting our 
data to Eq. (6), we establish that the incommensurate phase is 
critical and its central charge is c = 1 over the entire phase. 

While constructing the detailed phase diagram cross sec¬ 
tions, we found that while it was easy to approximate the 
locations of the phase boundaries, we often encountered dif¬ 
ficulties in precisely nailing down the central charge of the 
corresponding critical points. As an example, we note the ap¬ 


pearance of a few points with (apparently) high central charge, 
indicated by red color, on the direct topological-trivial phase 
boundary in Fig. 2(d). While in some cases there may be real 
physics associated to this behavior, we show in Appendix B 
that a primary source for these spurious effects is fitting to a 
region of the phase diagram that is just slightly off-criticality. 
We show that the central charge is very sensitive to the precise 
location of the critical point, and can easily give 0(1) errors 
even when only slightly tuned away from criticality, and even 
with reasonably large-size calculations. 

Additionally, although most phase boundaries were eas¬ 
ily identified, there are three regions where difficulties arise: 
(i) the trivial-incommensurate phase transition at = 0 and 
large 0 (lower-right corner of Fig. 2(d)), (ii) the topological- 
incommensurate phase transition at 0 = n/3 and small 6 
(upper-left corner of Fig. 2(e)), and (iii) the Lifshitz transition 
area for / = 0.5 and $ = Q ~ 7r/6as seen in Fig. 2(f). Re¬ 
gions (i) and (ii) are related by duality, and the explanation of 
the numerical difficulties in these regions may have a common 
origin. To explain, we recall that the trivial-incommensurate 
phase transition at 0 = 0 and large 0, i.e. region (i), is of the 
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(a) (b) (c) 

Figure 3. (Color online) Properties of the critical line at / = J — 1/2 for various values of 0 — f (a) Profile of the EE as a function of block 
size shows Lifshitz oscillations. We predict the oscillation length for f < 0.2 to be larger than our system size L = 200. (b) Energy gap 
between the ground and first excited state, which displays similar oscillatory behavior on varying the system size, (c) Characteristic oscillation 
lengths in the EE and energy gap, which are nearly identical for a large range of cf>. The (green) line is the fit with f = f~ 3 ' 75 + 1.16. 


Kosterlitz-Thouless type [42]. Hence, the correlation length 
decays as exp(c(T — Tkt ) -1 ^ 2 ) away from the transition 
point [61, 62], and this results in a long correlation length 
(compared to our system size L = 100) for this region of the 
phase diagram. The duality indicates that region (ii) may also 
be near a Kosterlitz-Thouless phase transition point. Thus, we 
attribute the issues with these regions as likely artifacts due to 
finite size effects. We elaborate further on this in Appendix C. 
The remaining region (iii) requires more discussion, to which 
we now turn. 

Lifshitz behavior — Let us now focus on the cross-section 
in Figs. 2(c), 2(f), which corresponds to f = 0. Since the 
system is self-dual on the line / = J, the trivial-topological 
phase boundary should just be the line / = J = 0.5, a fact 
verified in our numerical calculations when 6 = f are small. 
On top of the phase diagram we also plot the function / = 
[2 sin(30)] [1 + 2 sin(3</>)] -1 (in a thick circular dotted line), 
which represents an upper bound on the region in which exact 
parafermionic zero modes are expected to exist as proven in 
Ref. 13. The region of the phase diagram above this curve are 
guaranteed to not have exact parafermionic zero modes, de¬ 
spite still being in the topological phase with the topological 
ground state degeneracy. Along the critical line / = J = 0.5, 
c = 4/5 at the ferromagnetic point (</> = 0 = 0), and c = 1 
at the antiferromagnetic point = 0 = 7r/3) [50]. It is a pri¬ 
ori unclear how the central charge transitions from c = 4/5 
to c = 1, i.e., is it an abrupt jump at some transition point, 
or does it change incrementally in stages, or perhaps some¬ 
thing else entirely? Only a few studies address this question 
directly: among them is the work of Howes et al. [43] who 
used fermion analyses and series expansions to conjecture that 
a tricritical point connecting the ordered (topological), disor¬ 
dered (trivial), and incommensurate phases exists at exactly 
f = 0 = 7r/6. McCoy et al. [52, 53] studied the super in- 
tegrable line f = 0 = 7r/6 and suggested a modified picture 
with the incommensurate phase stretched all the way down to 
the point f = 0 = 0 and / = J = 0.5. Our results seem to 
support the latter picture, as we will further develop below. 


To address the questions posed above, we studied the crit¬ 
ical line f — J carefully. We observed (see Fig 3(a)) 
that before we reach the putative tricritcal (Lifshitz) point at 
f = 9 = 7r/6, the EE starts to show oscillatory behavior [63] . 
The frequency of the oscillations increases as we approach the 
Lifshitz point from small f = 6, and when further increasing 
f = 0 its amplitude dies out after the system clearly enters 
the incommensurate phase. Conventionally, a Lifshitz transi¬ 
tion point of this nature corresponds to a continuously vary¬ 
ing oscillation length, and in this case it is the length scale 
associated with the incommensurate order. Interestingly, the 
shapes of the EE oscillation curves match those observed re¬ 
cently in ID free, and interacting, fermion systems near Lif¬ 
shitz points where the Fermi surface is augmented by addi¬ 
tional Fermi points [44]. Thus, our result adds to the evidence 
of Ref. 44 that these types of EE oscillations are a finger¬ 
print of the Lifshitz-type phase transition. As an aside, we 
mention that the Lifshitz oscillations are only present in the 
EE when one uses open boundary conditions. One can easily 
check this by calculating the EE for free fermions as a func¬ 
tion of next-nearest neighbor hopping [44], but with periodic 
boundary conditions (shown in Appendix D). 

To quantitatively study the nature of this critical regime we 
want to investigate the variation of the central charge. How¬ 
ever, in the presence of oscillations in the EE, we must mod¬ 
ify Eq. (6) if we wish to extract the central charge. Empiri¬ 
cally, the observed oscillations appear to have a similar form 
to those in the work Ref. 64, and we propose a phenomeno¬ 
logical scaling form which can fit the EE with oscillations: 


L 


7 d 


S(1)cot = - In — Sin — + S 0 + 


cos(27rZ/C + p) 
(L/2 — \L/2 — l\) w ’ 


where the first two terms are the same as in Eq. (6), and the 
third term incorporates oscillations and a symmetrized damp¬ 
ing function. The parameter (is the oscillation length and p is 
a phase factor. These parameters, along with the exponent w, 
are free-parameters determined by fitting. Some representa¬ 
tive fits are shown in Figs. 4(a) and 4(b), which clearly capture 
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• • |_=100 c=0.8942 (=199.44 w=l.02 

• • |_=200 c=0.9070 (=131.21 w=0.66 

• • L=300 c=1.3929 (=192.06 w=0.67 

• . |_=400 C=1.4250 (=254.21 w=0.64 


c 0 200 400 

(a) 



(b) 




0 

(C) 


Figure 4. Panels (a) and (b) show the profile of the entanglement entropy (as a function of block size) for various values of system size at 
<j) — 6 — 0.25 and <j) — 0 — 0.35 respectively. Panel (c) shows the central charge obtained by fitting the entanglement entropy with the 
corrected formula along the line 4> = 6 and f — J — 0.5. The two dashed lines are at c = 0.8 and c— 1. The arrow indicates the trend of the 
peak when L is increased. 


the sub-leading oscillations accurately. 

The results of calculating the central charge from this pro¬ 
cedure are shown as a function of </> in Fig. 4(c). One can see 
that there is still an unaccounted for effect that leads to a peak 
in the central charge at a system-size dependent 4> value. More 
careful inspection reveals that the peak is located at a that 
corresponds to an oscillation length ( & L/2. Thus, as seen 
in the figure, the peak location </>*, occurs at values closer and 
closer to 0 = 6 = 0 when system size is increased, and all 
other parameters remain fixed. Our observations indicate that 
the central charge converges to c « 1 when </> > </>*, and 
c « 4/5 for (j) < (j)*. This strongly suggests that the transition 
from c = 4/5toc=l along the line f = J = 0.5 is an abrupt 
one that occurs at (j) = 6 = </>*. From our numerical data it 
appears that </>* — 0 as L oo. Hence, our data supports a 
scenario where there is an immediate onset of oscillations as 
one tunes away from <fr = 0 = 0 in the thermodynamic limit. 

We corroborate this by observing that oscillations are not 
seen in the EE if the oscillation length itself exceeds the sys¬ 
tem size L. For example, for L = 200, the oscillations are not 
explicitly visible for < 7r/12, however upon increasing the 
system size, with all other parameters fixed, the oscillations 
appear over a larger region of 0, as is shown in Fig. 4(a). As 
0 is decreased the oscillation length increases, and thus we 
must use larger and larger systems to observe the oscillations. 
Thus, we believe that this is evidence that, in the thermody¬ 
namic limit, the oscillations are a feature for all 0 = except 
0 = (j) = 0. An alternate scenario, which we can not rule 
out completely based on this numerical data, is that the in¬ 
commensurate phase persists to small but non-zero values of 
0 = (j). Thus, a conservative estimate of the location of the 
tricritical point isO < (0 = <p) < 0.25, which is well be¬ 
low the previously conjectured location at 0 = </> = 7r/6. We 
aim to shed further light on this transition through larger scale 
simulations in future work. 

Finally, we note that matching oscillations are observed in 
the splitting of the lowest two energy states (Fig. 3(b)), as 
a function of system size. We can extract the characteris¬ 


tic length scale ( of the oscillations from both the EE (for a 
given system length), and the energy gap (as a function of sys¬ 
tem length). Our results are shown in Fig. 3(c) where a clear 
correlation between the two is observed for 0 = 0 < 7r/4. 
The solid (green) line in Fig. 3(c) is the fit of the oscillation 
length for </> = 0 < 7r/4 to the function ( = <p~ 3 - 75 + 1.16. 
When = 0 = 0, the oscillation length appears to di¬ 
verge, indicating that no such oscillations survive in the non- 
chiral 3-state Potts model limit. Attempts to relax the fit with 
( = (0 — 0*) _7? + const (i.e. with a possibly non-zero (/>*) 
gave 0* ^ 0.09 indicating that the conjectured tricritical point 
may be in close proximity to =0. 

Conclusions — In summary, we have mapped out the three 
dimensional phase diagram of the Z 3 chiral clock model using 
the density matrix renormalization group method. Using the 
entanglement entropy (of the half-chain) as a diagnostic, we 
have been able to locate the phase boundaries of the various 
topological-trivial-incommensurate phase transitions. Quanti¬ 
tatively, we have also been able to see the variation of the cen¬ 
tral charge along the various critical surfaces that divide these 
phases. Another outcome of this study is the identification of 
the Fifshitz transition using the entanglement entropy, along 
with an estimate of the location of the putative tricritical point. 
We discussed several competing qualitative scenarios for the 
cross section of the phase diagram in which the tricritical point 
has been predicted to exist. Our data suggests that the tricriti¬ 
cal point (along / = J = 1/2) is not at (j) = 0 = 7r/6: rather 
we find it to be shifted to a much smaller value in the range 
0 < 0 = (j) < 0.25. 

Finally, our results must be viewed in a broader con¬ 
text as providing further confirmation of the stability of the 
parafermionic topological phase to chiral interactions, over a 
wide range of parameters. We expect a further study of this 
and related models to elucidate the conditions under which 
these phases can be practically realized. 
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These dual operators satisfy /i 3 = 1, u 3 = 1, and {iv = w vfi. 
The dual Hamiltonian is then 


L—l 


L 


^-dual 





—iO 


i0 + h.c. 

3 = 1 


(A9) 

Comparing with the original Hamiltonian, the dual Hamilto¬ 
nian returns to the original form if we exchange 0 and 0, and 
at the same time J and /. 


Appendix A: Properties of the 3-state chiral clock model 


The Hamiltonian in Eq. (1) has the same properties when 
either of the two phases 0 are shifted by multiples of |7r. 
To show this we can see that the transformation 


. 2rm . 2m7T 

O' -)> e + — </>' -> </> + — 

changes the Hamiltonian to: 

L L—l 


(Al) 


H 3 = -fu~ n Y, T]e-» - Jui~ m Y 
i = 1 j =i 

Then we can redefine the operators: 


+ie ld + h.c. 


(A2) 


r' =UJ n r a' 2j = uj rn a 2 j <T f 2 j+1 = CT2J+1- (A3) 

This new set of operators preserves the properties r 3 = a 3 = 
I, ar = u tcf , where cc = e 27 ™/ 3 . After this redefinition we 
end up with a Hamiltonian with the same form as the original. 

Additionally, the transformation that flips the signs of the 
two phases at the same time, i.e., 


O' ^-0 4>' -cj) (A4) 

changes the Hamiltonian to, 

L L -1 

Hs = -fY r 3 e ~^ - J X a i a Ui e ~ ie + h - c - ( A5 > 

3 = 1 3 =1 


Here we can redefine the operators as 


T — T ^ 


(j' = (T^ 


(A6) 


to recover the form of the original Hamiltonian. 

We can also just flip the sign of just one of the phases, say 
<// —>• — 0, and then the redefinition: 

T i = T -i a j = a -j ( A? ) 

leaves the Hamiltonian unchanged. If instead we flipped 
the sign of 0 , will need a transformation that involves both 
Eqs. (A6) and (A7). 

Finally, we can consider the duality transformation: 

3 

Mj+i = n rfc ’ v o+h =a j a j+ i ‘ ( a§ ) 

k=1 


Appendix B: Extracting Central Charge Near Critical Points 

When performing the fit to EE data obtained from a finite 
size system, and for a point in parameter space that is close 
to (but not at) a critical point, it is often difficult to obtain a 
reasonable estimate of the central charge. One possible expla¬ 
nation is that, when the system size is smaller than the corre¬ 
lation length, the fit to Eq. (6) may appear to be good, but the 
central charge obtained from the fit may not match the actual 
central charge of the nearby critical point. This is not unique 
to our model, and we were also able to observe this effect for 
free Dirac fermions with a tunable mass term as shown below. 
Eventually, if the system is tuned off criticality, and when the 
system size is larger than the correlation length, the EE will 
saturate and hence reveal the gapped phase. 

To provide an example of such behavior, we refer to known 
analytic results that the central charge should be 4/5 at (/ = 
J = 0.5, (j) = 0 = 0), and zero for all other / at <fi = 0 = 0. 
In Fig. 5, we show that at the critical point / = J = 0.5, the 
central charge is c = 0.81 ±0.01, close to the analytical re¬ 
sult. However, when we are slightly away from this point, say 
/ = 0.499, the system still appears critical with an (appar¬ 
ent) central charge of c = 1.58, much larger than the expected 
value of 0.80. On going slightly further away, / = 0.495, a 
plateau in the EE profile is seen consistent with our expecta¬ 
tion of a gapped phase. Thus, the fitting procedure produces 
misleading results in the neighborhood of the critical point, 
and can make it difficult to determine the central charge for 
critical points in which the position of the point is not known 
to extremely high accuracy. 


1. Near Critical behavior 

To further confirm our discussion above we performed sim¬ 
ilar calculations for ID gapless Dirac fermions using exact 
diagonalization. We use the 2-band free-fermion lattice Dirac 
model as the test model 

H = ~En ( ic Ui,t c n,i + * 4 + 1 ,^ + h - c *) 

— En ~ C n+l t 4 Cn >4 + h.G^j 

+ (2 — m) — 4,4^4) 

This model is gapless at k = 0 if m is zero, and the critical 
point should have a central charge of 1. If m is tuned away 
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l 


Figure 5. The EE as a function of the subsystem size l at 4> = 0 = 0 
and several different / close to or at the critical point (/ = 0.50). 
From the highest curve to the lowest one, the corresponding / is 
0.495, 0.499, 0.500, and 0.501. The central charges obtained from 
the fitting are shown in the legend. For / = 0.495 and f — 0.501, a 
plateau in the EE is seen indicating a gapped phase. For / = 0.499, 
an apparent critical phase is seen which is attributed to an artifact of 
finite size effects. 


from zero the system exhibits an energy gap of the size 2m. 
For our entanglement calculations the system was filled to half 
filling, such that when it is gapless, the filling hits exactly at 
the Dirac node, and if it is gapped, the filling includes all the 
states in the lower band. In this model the correlation length 
is controlled by the scale 1/m (with units it would b t hvp/m 
but h and vp are effectively unity for our model). 

To compare closely with our DMRG results we fit the cen¬ 
tral charge of this model using entanglement scaling with 
open boundary conditions. When gapless, we find the cen¬ 
tral charge to 2 or 3 digits of accuracy. For example, we 
find c= 1.006 when the chain is of length 400. In addition to 
calculating the scaling law over the entire chain we can im¬ 
prove the fit by taking symmetric cuts around the center of 
the chain which reduces the edge effects. We get slightly im¬ 
proved accuracy for ranges such as 120-280, i.e., c=1.004. If 
we increase system size to L=500 and fit over 120-380 we find 
c= 1.003. 



40-360 

120-280 

40-80 

L=300 

1.0198 

1.0202 

1.0199 

L=400 

1.022 

1.024 

1.021 

L=500 

1.0243 

1.0273 

1.0219. 


Table I. The central charge obtained by fitting from different region 
of the system and different system size L. The mass gap is set to be 
771=1/10000 

Now let us perturb the system slightly away from the crit¬ 
ical point. For this test we turn on a gap size of m= 1/10000 
as a start. As an estimate, this should give a correlation length 


of £ =10000 sites. For system size 400, if we fit from 40-360, 
we find c= 1.022; if we fit from 120-280 we find c= 1.024. If 
we try to fit a different range, e.g., 40-80 we find c=1.021. 
Either way, the result is already 1% different than the gap¬ 
less case even for this tiny gap (compared with the band¬ 
width). Next we repeated the same 3 fits for L=300 and we 
find c= 1.0198, 1.0202, 1.0199. And then for L=500 and find 
c= 1.0243,1.0273,1.0219. These results are summarized in Ta¬ 
ble. I. We observe that the fits get worse when we increase the 
system size, and when we fit over the region restricted mostly 
to lie over the center. The latter result may be expected since 
the scaling function varies most slowly over the center. The 
fact that the fits get worse as we increase system size is most 
likely just an indicator that there is a finite correlation length 
and that the critical scaling form will eventually break down. 
For additional tests we also fit the central charge for larger (but 
still very small) mass gaps with m = 1/1000 and m = 1/100 
in Table. II and Table. Ill respectively. 



40-360 

120-280 

40-80 

L =300 

1.1374 

1.1633 

1.1230 

L=400 

1.1699 

1.2082 

1.1372 

L =500 

1.2014 

1.2499 

1.1473 


Table II. The central charge obtained by fitting from different region 
of the system and different system size L. The mass gap is set to be 
771=1/1000 



40-360 

120-280 

40-80 

L =300 

1.8811 

1.7267 

1.9246 

L=400 

1.7563 

1.4033 

1.9438 

L =500 

1.5701 

1.0878 

1.9385 

L=600 

1.3874 

0.84507 

1.9273 


Table III. The central charge obtained by fitting from different region 
of the system and different system size L. The mass gap is set to be 
771=1/100 

We see that when we are tuned near, but not at, the criti¬ 
cal point the best fits in the gapped case seem to come from 
smaller system sizes, and over ranges which do not include the 
flat middle portion of the scaling range nor the far tails of the 
scaling range. The unfortunate thing is that once we are a bit 
further away from the critical point this optimized fitting pat¬ 
tern no longer works. In this case none of the fitting regimes 
we used give accurate results because the system begins to re¬ 
veal its gapped nature. We do find something close to c = 1 
when m = 1/100 and L = 500 (Table. Ill), but this seems ac¬ 
cidental since we tested it for L=600 and got a worse results. 
From this data we would claim that for the Dirac model when 
the central charge differs by 20% from its expected value then 
we are too far away from the critical point to do any fitting and 
should claim that it is not critical. In fact for a system size of 
500 and mass gap of mm 1/1000 the fitted values are closer 
to 6/5 instead of 1 and could easily lead to misidentification 
of critical points in models where their location is not known 
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exactly. 



x 

Figure 6. The transformed curve of entanglement entropy. Here x is 
defined as ln( ^ sin ^), where l is the block size and the system size 
L & 200. For small x the curves are ordered by increasing the mass 
gap (lowest mass is the lowest curve). 

As a possible diagnostic we plot the entanglement entropy 
as a function of x = ln(^ sin ^), where l is the sub-block 
size. The slope of the entanglement entropy vs x should be 
interpreted as c/6. The only feature that could be used as 
a diagnostic is that if the transformed curve has a decrease 
in slope then we are definitely too far away from a critical 
point to fit properly as can be seen in Fig. 6. The final two 
curves have clear decreases in slope as we move far away from 
criticality. Note that all these artificially high central charges 
only occur when we use open boundary conditions. As can 
be seen in Fig. 7, the central charge first goes up then drops 
for open boundary conditions when we tune the system away 
from criticality. However, it decreases monotonically for pe¬ 
riodic boundary conditions. 



Figure 7. The central charge obtained by fitting the entanglement 
entropy from site 40-160 for a 200 sites chain. The blue open circles 
are for open boundary conditions and the red dots are for periodic 
boundary conditions. 


Appendix C: Kosterlitz Thouless transition in the 3-state chiral 
clock model phase diagram 

In the main text, we studied several 2D cross sections of 
the 3D (/, J = 1 — /, 0, (j)) phase diagram of the chiral clock 
model. The 2D cross sections corresponding to 0 = 0 (see 
Figs. 2(a) and 2(d)) and 0 = 7t/3 (see Figs. 2(b) and 2(e)) 
showed some regions whose phase boundaries could not be 
located. This was attributed to finite size errors, which we 
now address. 



/ 


(a) 



/ 


(b) 

Figure 8. Panel (a) shows the entanglement entropy (for the central 
cut), as a function of / for 6 = 1.00. The EE increases for larger 
sizes of the system for / from 0.07 to 0.17, indicating a critical phase 
at this region, (b) shows the corresponding central charge calculated 
for various system sizes. The change of the central charge becomes 
sharper for larger systems. 

We first discuss the features seen in Figs. 2(a) and 2(d), i.e., 
the cross section for 0 = 0. For small / and large 0, the phase 
transition between the topological and trivial phase is indirect: 
it is mediated by the incommensurate phase. To establish the 
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(a) 



(b) 

Figure 9. The profile of the entanglement entropy as s function of 
block size l at 4> = 7t/3, 9 = 0 and (a) / = 0.8 (b) / = 0.9 for 
different system size. The continuous lines are the fit to the DMRG 
data. 


fact that the incommensurate region is of non-zero extent, we 
performed finite size analyses on both the entanglement en¬ 
tropy and central charge as is shown in Figs 8(a) and 8(b) This 
extent is found to be from / ^ 0.07 to / ~0.15. We find that 
the central charge of the trivial-incommensurate transition is 
consistent with that of the Kosterlitz-Thouless (KT) type, i.e., 
c = 1 [42]. 

Because of the duality in the Hamiltonian (Eq. 1), the phase 
diagram is symmetric with respect to the line / = J = 0.5, 
(j) = 6. Thus, the above mentioned phase transition is dual to 
the incommensurate-topological phase transition, for large 0 
and small 0. That is to say, the region with the smooth change 
of the central charge in the lower-right corner of Fig. 2(d) is 
dual to the (red) region in the upper-left corner of Fig. 2(e). 
This region, being near the KT phase transition point is also 
plagued by finite size errors: the correlation length is long 
compared with the system size (L = 100). 

To test this assertion, we studied the (apparently) large cen¬ 
tral charge that was calculated near the critical region, as is 
shown in Fig. 9. For example, as is shown in Fig. 9(a), the 
point (j> = 7 t / 3 , 6 = 0, and / = 0.8 appears to be critical, but 
for larger system sizes is shown to be gapped. We base this 
conclusion on the appearance of a saturation plateau in the 


profile of the EE scaling as a function of subsystem size. As 
a comparative check, we went deeper into the critical regime 
(i.e. / = 0.9). As can be seen in Fig. 9(b) and as is expected, 
we found no such plateau in the EE. 


Appendix D: Lifshitz Transition for free fermions 


For comparison with our discussion of the Fifshitz transi¬ 
tion in the chiral clock model we consider a version with ID 
free fermions hopping on a chain with nearest neighbor and 
next nearest neighbor hopping. As the n.n.n. hopping is in¬ 
creased additional Fermi-points can enter the spectrum and 
eventually hit the chemical potential which leads to a Fifshitz 
transition of the Fermi-surface topology. As our model we 
consider free fermions with next nearest neighbor hopping. 


H = - 




[4+l C n + tc l+ 2 c n + h.C. 


(Dl) 


Here, t is the parameter for the next nearest hopping. The 
energy spectrum of this model is E = — 2 cos(fc) — 2t cos(2 k). 
When t increases from zero, the topology of the Fermi surface 
at zero energy changes from two points to four points at t = 1, 
which is the Fifshitz transition. 

We calculate the entanglement entropy of this model with 
open boundaries and the periodic boundaries at half filling. 
The results are shown in Fig. 10 and Fig. 11 and one can im¬ 
mediately recognize the pattern of oscillations that we saw 
earlier for the chiral clock model. One interesting thing to 
notice is that the oscillations go away when we use periodic 
boundary conditions. This model, and the related entangle¬ 
ment properties, are carefully studied in Ref. 44. For periodic 
boundary conditions the curves gradually increase from the 
scaling form with c = 1 to a scaling form of c = 2 which is 
the result expected for two sets of left and right movers at the 
Fermi-level. 
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